NASA Technical Memorandum 100089 


Potential Applications of 
Computational Fluid Dynamics 
to Biofluid Analysis 

D. Kwak, Ames Research Center, Moffett Field, California 
J. L. C. Chang, Rocketdyne, Rockwell International, Canoga Park, California 
S. E. Rogers, Sterling Federal Software, Palo Alto, California 
M. Rosenfeld, Ames Research Center, Moffett Field, California 


April 1988 


NASA 

National Aeronautics and 
Space Administration 

Ames Research Center 

Moffett Reid, California 94035 




POTENTIAL APPLICATIONS OF COMPUTATIONAL 
FLUID DYNAMICS TO BIOFLUID ANALYSIS 


D. Kwak,* J. L. C. Chang, t S. E. Rogers,}: and M. Rosenfeld § 
NASA Ames Research Center, Moffett Field, CA 94035, U.S.A. 


SUMMARY 

Computational fluid dynamics has been developed to the stage where it has become an 
indispensable part of aerospace research and design. In view of advances made in aerospace 
applications, the computational approach can be used for biofluid mechanics research. In 
the present paper, several flow simulation methods developed for aerospace problems are 
briefly discussed for potential applications to biofluids, especially to blood flow analysis. 

INTRODUCTION 

To date major advances in computational fluid dynamics (CFD) have been made in 
aerospace engineering. With the advent of supercomputers as well as the development of 
fast algorithms, computational flow simulations have become a practical means for 
aerospace designs. Therefore, it will be of considerable benefit to medical researchers to 
extend this CFD technology to biofluid analysis. One of the important areas of biofluid me- 
chanics deals with blood flow, such as heart and blood vessel problems, which is relevant to 
cardiovascular diseases and their treatment. Understanding the flow phenomena by numer- 
ical simulations will be of significant value toward finding treatments for these problems. 
Therefore, the potential payoff to human health will be tremendous. 

Blood flow is very complicated in many respects: The fluid may exhibit significant 
non-Newtonian characteristics locally and the geometry is usually very complicated. For 
instance, the human aorta has large curvatures combined with very irregular lumen cross 
sections and the walls are elastic and change shape (about a 10% increase in diameter during 
the systole). In an artificial organ, as red blood cells go through high-shear turbulence 
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regions, they may be damaged; the downstream region of an artificial heart valve is an 
example. The flow is unsteady, possibly periodic, and very viscous and incompressible. 
This problem is very much interdisciplinary and an attempt for a complete simulation would 
be a very formidable task. However, an analysis based on a simplified model may provide 
much needed physical insight into the blood flow analysis. For a more comprehensive study 
on blood flow, a large number of publications are available (refs. 1-10) and is not reviewed 
here. 


In the recent past, viscous incompressible flow solvers have been developed at NASA 
Ames Research Center. This research was motivated from realistic needs for 
three-dimensional simulations of aerospace applications (refs. 11-14). Naturally, compu- 
tational efficiency has been of primary importance in addition to accuracy and robustness. 
The formulation is based on a Newtonian fluid assumption. However, since the govern- 
ing equations are solved in a generalized coordinate system, viscosity is allowed to vary in 
space and time. These flow solvers can be applied to current blood flow problems, such 
as flow through an artificial heart (ref. 9), pulsatile flow in arterial bifurcations (ref. 6) and 
flow in aneurysms (ref. 7). A full simulation of viscoelastic flow is very difficult because of 
the nonlinearities of the fluid (ref. 5). However, as a first step toward full simulations, non- 
Newtonian effects of the blood flow can be simplified by a constitutive model for viscous 
stresses. 

In the present paper, viscous incompressible flow CFD methods will be discussed for a 
potential extension to blood flow problems. The focus is on flow solvers recently developed 
by the authors, and on relevant features pertaining to blood flow simulations. 

FORMULATION 


Unsteady, three-dimensional, viscous, incompressible flow with constant density is 
governed by the following Navier-Stokes equations: 
continuity equation 



( 1 ) 


momentum equation 


duj | dutU; _ dp [ dr y 
dt dxj dxi dxj 

where, t is the time, i, the Cartesian coordinates, Uj the corresponding velocity components, 
p the pressure, and Ty the viscous stress tensor. Here, all variables are nondimensionalized 
by a reference velocity and length scale. 
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To close the governing equations (1) and (2), the viscous stress tensor is modeled by a 
constitutive relation. For a Newtonian fluid, this can be written as 


T\j — 2 vS%j Rxj 


(3) 


Here, v is the kinematic viscosity, Sij is the strain rate tensor, and R,j is the Reynolds 
stresses. The strain rate tensor is defined by 



duj 

dxj 


(jX i 


(4) 


Various levels of closure models for Rij are possible (ref. 15). If turbulence is simulated by 
an eddy viscosity model, the viscous stress tensor can be replaced by the following form: 


T\j = 2 vrSij (5) 

where ut is the effective total viscosity. Therefore, the present formulation allows vari- 
able viscosity which can be represented by a time-dependent constitutive relation. Non- 
Newtonian effects of blood can be included in a simple manner by an equation of state for 
vt- In a more sophisticated approach, a currently used differential model (ref. 5) for r i} ; 
namely, 

= f(t,p, Sij, material) (6) 

can be incorporated in the analysis. This equation can be solved in a decoupled mode from 
the governing equations by lagging the time advancement by one step. In the remainder of 
this paper, our discussion will be limited to an algebraic expression as in Equation (5). 

The geometric variation for blood flow computations is diverse and naturally the com- 
putational approach has to take this into account. To perform calculations on 
three-dimensional, arbitrarily shaped geometries, the following generalized independent 
variables are introduced which transform the physical coordinates into general curvilinear 
coordinates 
r = t 

, * (?) 
« = c( x,y, z,t) 

Applying the transformation to the governing equations yields a generalized form of the 
incompressible Navier-Stokes equations cast in curvilinear coordinates. Various simplifi- 
cations can be made on these equations depending on the geometry and the boundary-layer 
thickness. For instance, in many problems of external aerodynamics where the viscous 
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region is confined to a thin region, the Navier-Stokes equations are often simplified by a 
thin-layer approximation or by a parabolic form. For internal flows, such as the blood flow 
in the human body, viscous effects are important in the entire flow field. Therefore, the set 
of governing equations need to be solved in their complete form. 

SOLUTION METHODS 

A major problem in solving the incompressible Navier-Stokes equations comes from 
the lack of a pressure term in the continuity equation. In realistic three-dimensional appli- 
cations, satisfying continuity in a reasonable amount of computing time becomes a primary 
issue. Various methods of resolving this problem have been developed (ref. 15). However, 
the present discussion is limited to a primitive variable formulation in generalized curvi- 
linear coordinates. In this section, three flow solvers recently developed by the authors are 
summarized to suggest a way of extending this type of computational tool to blood flow 
simulations. Derivation of equations and algorithmic details can be found in full-length 
versions of the authors’ publications (refs. 11-14). Other numerical methods developed 
in the past for simulating viscous incompressible flows can be found in the literature (see 
(ref. 15) for a review). 


Pseudocompressibility method 


Recent advances in the state of the art in CFD have been made in conjunction with com- 
pressible flow computations. Therefore, it is of significant interest to be able to use some 
of these compressible flow algorithms. To do this, the artificial compressibility method 
(ref. 16) can be used. In this formulation, the continuity equation is modified by adding a 
time derivative of the pressure term resulting in 


1 dp dui 
fi dt dxi 


( 8 ) 


where fi is a pseudocompressibility coefficient. Together with the unsteady momentum 
equations, this forms a hyperbolic-parabolic type of time-dependent system of equations. 
This method was originally intended for the steady-state computation of incompressible 
flow (refs. 11 and 12). However, by introducing a pseudo-time iteration, this can be made 
time-accurate (refs. 14 and 17). 


By constructing a pseudocompressible form of the governing equations, 
fast, implicit schemes developed for compressible flows, such as the 
approximate-factorization scheme (ref. 18), can be implemented. Based on this approach. 
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a computer code (INS3D (ref. 11)) has been written for obtaining steady-state solutions. 
The spatial discretization uses second-order central differencing with additional numerical 
dissipation terms. This code has been validated and applied to numerous three-dimensional 
problems (refs. 11, 12, 19, and 20). 

To obtain time-accurate solutions using the pseudocompressibility formulation, it is 
necessary to satisfy continuity at each time step by subiteration in pseudo-time. To use 
a large time step in the pseudo-time iteration, an upwind differencing scheme based on 
flux-difference splitting is used combined with an implicit line relaxation scheme. This 
removes the factorization error and the need for numerical dissipation terms. A second 
code (INS3D-UP) based on this method has been validated and excellent results have been 
obtained (ref. 14). 


Fractional-step method 

The fractional-step method is another approach that can be used for time-dependent 
computations of the incompressible Navier-Stokes equations. Here, the discretized equa- 
tions are advanced in time by decoupling the solution of the momentum equation from that 
of the continuity equation. The common application of this method is done in two steps. 
The first step is to solve for an auxiliary velocity field using the momentum equation in 
which the pressure-gradient term is approximated by its value at the previous time step. 
In the second step, the pressure, which maps the auxiliary velocity onto a divergence-free 
velocity field, is computed. Various other operator splittings can be adopted by treating the 
momentum equation as a combination of convection, pressure, and viscous terms. 

A third flow solver (INS3D-FS) based on this approach using a staggered grid has been 
developed (ref. 13). The governing equations are discretized conservatively using finite 
volumes. Rather than choosing the Cartesian velocity components as dependent variables, 
the volume fluxes over the faces of the computational cells are used. They are equivalent 
to the contravariant velocity components described in a staggered grid. This procedure, 
combined with accurate and consistent approximations of the geometric quantities, satisfies 
the discretized mass conservation equation. In the second step, a novel four-color scheme 
is devised for solving the resulting Poisson equation for the pressure correction. Several 
computational results have been compared with experiments and other numerical solutions 
in reference 13. 
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COMPUTATIONAL PROCEDURE 


Grid generation 

Commonly used grid topologies can be classified basically into three different types, 
namely, O-, C- and H-grids. In blood flow applications, the blood vessels have tubular 
cross sections and can branch off into several smaller tubes. The flow can go through ir- 
regular passages and may encounter rapid turns and expansions as in an artificial heart 
(refs. 8 and 9). Naturally, computational simulation of these flows needs to be performed 
in more than one zone that are separated by geometric characteristics. For a multiple zone 
calculation, geometric continity is desirable for avoiding complicated interface procedures. 
Various grid generators are available in the literature. 

To illustrate a possible way of generating a realistic grid, the geometry of a glass 
aneurysm model, proposed by Liepsch et al. (ref. 7), is chosen. Note that, by combin- 
ing the two H-grids, the grid in the cavity is continuous with the grid in the main blood 
vessel. Figure 1(a) shows that the grid for the physical domain and figure 1(b) shows the 
grid for the computational domain (generated by L. Chang: unpublished note). Here, the 
circular cross sections are mapped by an H-grid. This topology is very flexible in cluster- 
ing the grid points in the area of interest; however, it does so at the expense of introducing 
singular points. These singularities are hidden in the corner of the solid boundary and can 
easily be handled by explicit boundary conditions. 

Computed examples 

The three flow solvers just described have been validated by computing basic fluid dy- 
namics problems (refs. 11-13). A few examples are listed in this section for demonstrating 
the capability of those solvers for potential applications to blood flow problems. 

Bifurcation is an important blood flow problem, since a low-speed recirculating region 
and local regions of high shear stress may play an important role in the formation of atherotic 
plaques and thrombi as well as releasing hemoglobin in the blood stream. As an idealization 
of a branching human circulatory system, a numerical solution of a 90° bifurcation problem 
(ref. 6) is shown in figure 2. This result is obtained using the INS3D-UP code (ref. 15). A 
more realistic geometry can be modeled by a grid similar to the aneurysm shown in figure 1 . 
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In figure 3, the time evolution of a two-dimensional cavity flow is shown by time 
dependent stream function contours. This result is obtained by the INS3D-FS code (ref. 1 3). 
This example illustrates the potential of extending the flow solver to an aneurysm analysis. 

As an example of how the present CFD capabilities are being used in real world ap- 
plications, work related to the Space Shuttle main engine (SSME) flow simulation (refs. 19 
and 20) is briefly discussed here. The geometric complexity of the SSME power head is 
comparable to blood flow problems, such as flow in an artificial heart. In the SSME staged 
combustion cycle, the fuel is partially burned at very high pressure and relatively high tem- 
perature in the prebumers. The resulting hot gas is used to drive the turbine and is then 
discharged from the gas turbine to the annular turnaround duct with a 180° U-turn before it 
diffuses into the fuel bowl. Flow through this assembly was simulated using INS3D code. 
The flow in the U-tum is especially interesting when the flow is turbulent because the tur- 
bulence structure is greatly influenced by streamline curvature. Therefore, in this case, it 
is necessary to incorporate strong curvature effects in the turbulence modeling. For the 
present paper, a length scale determined by the point of minimum vorticity is incorporated 
into an extended Prandtl-Karmann mixing-length theory. The combination of these auto- 
matically account for the curvature effect. Full details of this model are given in reference 
(ref. 21). In figure 4(a), the grid in the 180° turn region is shown. In figures 4(b)-4(i), com- 
puted velocity profiles using this algebraic turbulence model at Re=10 5 are compared with 
experiment (ref. 22). Considering the difficulties associated with curvature effect modeling, 
the results are very promising. 


Postprocessing 

Three-dimensional simulation of realistic flow produces an avalanche of data. There- 
fore, fast postprocessing tools are necessity for analyzing the data. Various 
three-dimensional graphic softwares have been developed in parallel with supercomputers. 
Graphics work can be performed on a work station independently or interactively with the 
main supercomputer. This provides an invaluable means in interpreting and understanding 
the results of the computation. Further details of the recent development in postprocessing 
can be found in the references cited (refs. 23 and 24). 

CONCLUDING REMARKS 

This paper presents a potential procedure of applying incompressible Navier-Stokes 
flow-solvers to blood flow problems. Computed examples on several generic problems 
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relevant to blood flow are presented. For problems related to large blood vessels, the present 
Newtonian formulation is an acceptable approximation for gaining insight into the flow 
physics involved. However, for a more complete and accurate analysis, it will be required 
to model the non-Newtonian effect as well as the transition and turbulence effects. Despite 
the limitation on physical models, the successful application of the incompressible Navier- 
Stokes solvers to the SSME analysis and redesign, provides an example of present CFD 
capabilities. These capabilities can be integrated into the blood flow analysis, since similar 
approximations to the physics can be made in that case. 
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Fig. 2 Numerical solution of 90° bifurcation (Re=496, vertical flow = 44 % of total flow): 
a) velocity vectors, b) stream function contours. 
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Fig. 3 Time evolution of 

at Re=1000. 
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